disp('lfG');

% t0{1} = 4;
% t0{2} = [2, -1];
% t = fn_simplify_oneterm(t0);
% 
% disp(t0{1});
% disp(t0{2});
% disp(t{1});
% disp(t{2});
% disp('hqlE');
% disp(fn_cal_one([2, -2, 2, 2]));
% t = fn_cal_one([2,2]);
% 
% one_term = {1,[2,3]};
% 
% T = {{1,[1,-2]},{2,[2,1,-1,1]},{4,[-4,4,-2,1]}, ...
%     {3,[-4,4,-4,4,-2,1]}};
% 
% T2 = fn_simplify(T);
% 
% T3 = fn_combine(T2);
% 
% fn_print(T);
% fn_print(T2);
% fn_print(T3);
% 
% T4 = fn_add(T3, T3);
% fn_print(T4);
% 
% T5 = fn_multiply(-0.15, T4);
% fn_print(T5);
% 
% t1 = {1.5, [1]};
% t2 = {-2, []};
% T6{1} = fn_prod_oneterm(t1, t2);
% fn_print(fn_simplify(T6));
% 
% T7 = fn_prod(T6, T6);
% fn_print(T7);
% 
T1 = {{1,[-1]},{2, [-2,2]},{3,[]}};
T2 = {{3,[-1,1]},{4, [2,-2]},{5,[]}};
fn_print(T1);
fn_print(T2);
T12 = fn_prod(T2, T1);
fn_print(T12);

T123 = fn_normal_order(T12, []);
fn_print(T123);

Num = 20;
Bqd = cell(1, Num-1);
Bq = cell(1, Num-1);

for q = 1:(Num-1)
    
   Bqd{q} = fn_generate_Bqd(q, Num);
   Bq{q} = fn_conjugate(Bqd{q});
    
end


fn_print(Bqd{1});
fn_print(Bq{1});
fn_print(fn_commutation(Bq{2}, Bqd{2}));
fn_print(fn_commutation(Bq{1}, Bqd{2}));
fn_print(fn_commutation(Bq{2}, Bqd{1}));
fn_print(fn_commutation(Bq{5}, Bqd{3}));
fn_print(fn_commutation(Bq{1}, Bqd{4}));

ck = cell(1, Num);
for q = 1:Num

    ck{q} = {{1, -q}};
    
end

fn_print(fn_commutation(ck{5}, Bqd{7}));

% T = fn_commutation(Bq{1}, Bqd{2});
% fn_print(T);
% T = fn_commutation(Bq{2}, Bqd{1});
% fn_print(T);
% 
% T = fn_commutation(Bq{3}, Bq{4});
% fn_print(T);
% T = fn_commutation(Bqd{1}, Bqd{2});
% fn_print(T);
% 
% T = fn_commutation(Bqd{2}, Bq{2});
% fn_print(T);

x = rand();
pshix = {{0,[]}};
for k = 1:Num
    
   pshix = fn_add(pshix, {{exp(-1j*k*x),-k}}); 
    
end

fn_print(pshix);

p = 3;
T = fn_commutation(pshix, Bqd{p});
T = fn_minus(T,fn_multiply(exp(-1j*p*x), pshix));
T = fn_simplify(T);
fn_print(T);

fn_print(fn_prod(T, Bqd{2}));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function T = fn_commutation(T1, T2)

T = fn_minus(fn_prod(T1, T2), fn_prod(T2, T1));
T = fn_normal_order(T, []);

end

function T = fn_minus(T1, T2)

T = fn_add(T1, fn_multiply(-1, T2));

end


function Bqd = fn_generate_Bqd(q, Num)

Bqd = {{0,[]}};

for k = 1:(Num-q)
    
   Tadd = {{1, [q+k, -k]}};
   Bqd = fn_add(Bqd, Tadd);
    
end

end


function t = fn_conjugate_oneterm(t)

a = t{1};
vx = t{2};

a = conj(a);
vx = -vx(end:-1:1);

t{1} = a;
t{2} = vx;

end


function T = fn_conjugate(T)

Num = length(T);

for n = 1:Num
    
    T{n} = fn_conjugate_oneterm(T{n});
    
end

end


function T = fn_normal_order(T0, v_normal_index)

Num = length(T0);

T = {{0,[]}};

for n = 1:Num
    t = T0{n};
    Tadd = fn_normal_order_oneterm(t, v_normal_index);
    T = fn_add(T, Tadd);
end

end


function T = fn_normal_order_oneterm(t, v_normal_index)

a = t{1};
vx = t{2};

if isempty(vx)
    T{1} = t;
    return;
end

T = {{1, []}};
while ~isempty(vx)
    
    vt = vx(abs(vx) == abs(vx(1)));
    
    if length(vt) == 2       
        if (vt(1) < 0 && ismember(abs(vt(1)), v_normal_index)) ...
                || (vt(1) > 0 && ~ismember(abs(vt(1)), v_normal_index))
            vt = vt(2:-1:1);
            T2 = {{1,[]},{-1,vt}};
        else          
            T2 = {{1, vt}};         
        end    
    else
        T2 = {{1, vt}};
    end
    
    T = fn_prod(T, T2);
    vx(abs(vx) == abs(vx(1))) = [];
    
end

T = fn_multiply(a, T);

end


function T = fn_prod(T1, T2)

N1 = length(T1);
N2 = length(T2);

T = cell(1, N1*N2);

n = 0;

for n1 = 1:N1
    for n2 = 1:N2
        t1 = T1{n1};
        t2 = T2{n2};
        t = fn_prod_oneterm(t1, t2);
        
        n = n + 1;
        T{n} = t;

    end

end

T = fn_simplify(T);

end


function t = fn_prod_oneterm(t1,t2)

a1 = t1{1};
a2 = t2{1};
v1 = t1{2};
v2 = t2{2};

v = [v1, v2];
a = a1*a2;

t{1} = a;
t{2} = v;

t = fn_simplify_oneterm(t);

end


function T = fn_multiply(a, T)

for n = 1:length(T)
    
   T{n}{1} = a*T{n}{1}; 
    
end

T = fn_simplify(T);

end


function T = fn_add(T1, T2)

T = [fn_simplify(T1), fn_simplify(T2)];
T = fn_simplify(fn_combine(T));

end


function res = fn_is_equal_term(t1, t2)

vx1 = t1{2};
vx2 = t2{2};

if length(vx1) == length(vx2)
    if sum(abs(vx1 - vx2)) == 0
        res = true;
        return;
    end
end

res = false;

end


function T = fn_combine(T)

ind = 1;

while ind<=length(T)
    
    vequal = zeros(1,length(T));
    vequal(ind) = 1;
    for p = (ind+1):length(T)
        if fn_is_equal_term(T{ind}, T{p})
            vequal(p) = 1;
        end
    end
    
    amp = 0;
    ind_non_zero = find(vequal == 1);
    for i = 1:length(ind_non_zero)
        t = T{ind_non_zero(i)};   
        amp = amp + t{1};     
    end
    
    t = T{ind};
    t{1} = amp;
    T{ind} = t;
   
    T(ind_non_zero(2:end)) = [];
    
    ind = ind + 1;
    
end

T = fn_kill_empty(T);


end



function T = fn_kill_empty(T0)

Num = length(T0);

va = zeros(1, Num);
for n = 1:Num

    one_term = T0{n};
    va(n) = one_term{1};
    
end

T = T0;
T(abs(va) < 1e-12) = [];
     
end


function fn_print(T)


fprintf(' result = \n');

if isempty(T)
    fprintf('0 \n');
end

Num = length(T);

for n = 1:Num
    
   t = T{n};
   fn_print_oneterm(t);
   fprintf('\n');
    
end

fprintf('===================== \n');

end


function fn_print_oneterm(t)


a = t{1};
vx = t{2};

if a>0
    fprintf('+');
end

fprintf('%4.2f * ', a);

if isempty(vx)
    fprintf('Id');
end

for i = 1:length(vx)
    
    n = vx(i);
    if n>0
        fprintf('C^%d ', n);
    else
        fprintf('C%d ', -n);
    end
    
end

end


function T = fn_simplify(T)

Num = length(T);

for n = 1:Num
    one_term = T{n};
    one_term_simplify = fn_simplify_oneterm(one_term);
    T{n} = one_term_simplify;
end

T = fn_combine(fn_kill_empty(T));

end

function t = fn_simplify_oneterm2(t)

a = t{1};
vx = t{2};

p2 = length(vx);

if p2 == 0
    return;
end

% disp(vx);
while 1
    
    n = vx(p2);
    p1 = find(abs(vx) == abs(n), 1);
    
    vx0 = vx(p1:p2);
    vx0 = fn_cal_one(vx0);
    
    if isempty(vx0)
        t = {0, []};
        return;
    end
    
    vx = [vx(1:(p1-1)), vx0, vx(p2+1:end)];
    
%     disp(vx);
    p2 = p1 - 1;
    
    if p2 == 0
        break;
    end
    
end

t{2} = vx;

end



function res = fn_simplify_oneterm(t)

vx = t{2};

Num = length(vx);

phase = 0;
for i = 1:Num
    for j = 1:(Num-i)
        
        x1 = vx(j);
        x2 = vx(j+1);
        
        if abs(x1)>abs(x2)
            
            vx(j) = x2;
            vx(j+1) = x1;
            
            phase = phase + 1;
            
        end
        
    end
end

res{1} = t{1}*(-1)^phase;
res{2} = vx;

res = fn_simplify_oneterm2(res);

% if isempty(res{2})
%     res{1} = 0;
% end

end


function res = fn_prod_single(t1, t2)

a1 = t1{1};
a2 = t2{1};

v1 = t1{2};
v2 = t2{2};

vx = [v1, v2];

t{1} = a1*a2;
t{2} = vx;

t = fn_simplify_oneterm(t);

end




function res = fn_cal_one(vx)

Num = length(vx);

for i = 1:(Num-1)
   
    x1 = vx(i);
    x2 = vx(i + 1);
    
    if x1 == x2
        
        res = [];
        return;
       
    end
    
end

if length(vx) <=2
    res = vx;
else
   if mod(length(vx), 2) == 0
       res = [vx(1), vx(end)];
   else
      res = vx(1); 
   end
end
    
    
end
